Application of QUBO solver using black-box optimization to structural design for resonance avoidance

Quadratic unconstrained binary optimization (QUBO) solvers can be applied to design an optimal structure to avoid resonance. QUBO algorithms that work on a classical or quantum device have succeeded in some industrial applications. However, their applications are still limited due to the difficulty of transforming from the original optimization problem to QUBO. Recently, black-box optimization (BBO) methods have been proposed to tackle this issue using a machine learning technique and a Bayesian treatment for combinatorial optimization. We propose a BBO method based on factorization machine to design a printed circuit board for resonance avoidance. This design problem is formulated to maximize natural frequency and simultaneously minimize the number of mounting points. The natural frequency, which is the bottleneck for the QUBO formulation, is approximated to a quadratic model in the BBO method. For the efficient approximation around the optimum solution, in the proposed method, we probabilistically generate the neighbors of the optimized solution of the current model and update the model. We demonstrated that the proposed method can find the optimum mounting point positions in shorter calculation time and higher success probability of finding the optimal solution than a conventional BBO method. Our results can open up QUBO solvers’ potential for other applications in structural designs.


Application of QUBO solver using black-box optimization to structural design for resonance avoidance Tadayoshi Matsumori * , Masato Taki & Tadashi Kadowaki
Quadratic unconstrained binary optimization (QUBO) solvers can be applied to design an optimal structure to avoid resonance. QUBO algorithms that work on a classical or quantum device have succeeded in some industrial applications. However, their applications are still limited due to the difficulty of transforming from the original optimization problem to QUBO. Recently, black-box optimization (BBO) methods have been proposed to tackle this issue using a machine learning technique and a Bayesian treatment for combinatorial optimization. We propose a BBO method based on factorization machine to design a printed circuit board for resonance avoidance. This design problem is formulated to maximize natural frequency and simultaneously minimize the number of mounting points. The natural frequency, which is the bottleneck for the QUBO formulation, is approximated to a quadratic model in the BBO method. For the efficient approximation around the optimum solution, in the proposed method, we probabilistically generate the neighbors of the optimized solution of the current model and update the model. We demonstrated that the proposed method can find the optimum mounting point positions in shorter calculation time and higher success probability of finding the optimal solution than a conventional BBO method. Our results can open up QUBO solvers' potential for other applications in structural designs.
Computing algorithms and hardware that aim to solve a quadratic unconstrained binary optimization (QUBO) have been recently developed. Quantum annealing (QA) 1,2 is one of the heuristic optimization algorithms for QUBO. QA utilizes quantum physics 3 to search for the optimal solution. QA is in the spotlight; after implementing QA in the quantum computer 4 , it was applied to some industry applications 5 . In the last few years, heuristic optimization algorithms that work on a classical computer, e.g., FPGA and GPU, have also been proposed for QUBO [6][7][8][9][10][11] . At present, these algorithms deal with a large number of design variables compared to QA that works on the quantum computer, and, in such a case, overcomes QA in terms of the time to solutions 12 . In the present paper, we call for an algorithm to solve QUBO, a "QUBO solver" without distinguishing whether the algorithm performs on a classical or quantum device.
QUBO solvers succeeded in solving some combinatorial optimization problems which are not directly formulated in a QUBO. It is known that some combinatorial optimization problems can be transformed into QUBO 13 . In addition, other transformation techniques for QUBO formulation have also been proposed, such as encoding techniques to a binary variable 14,15 and converting methods from higher-order polynomial to a quadratic model 16 . And, furthermore, many applications have been reported as the success cases of the QUBO formulation, including vehicle path generation 12,17 , traffic signal control 18 , portfolio optimization 19 , quantum chemistry 20 , and machine learning 21 .
When a QUBO solver is applied to a structural design in the industry, the QUBO formulation becomes a bottleneck. Simulations of differential equations based on the finite element method are widely used to design products accounting for shapes, materials, geometric constraints, and so on. When the design problem is formulated as the optimization problem, the differential equations and their solutions are dealt with by constraint and objective functions, respectively. For example, in the topology optimization which is a structural design method 22,23 , design variables represent the structure shape, and a partial differential equation constrained optimization is formulated. The bottleneck of the QUBO solver in the design optimization is that the constraints and objectives cannot be explicitly expressed as the quadratic polynomial of binary design variables. Design methods  [24][25][26][27] . The BBO encodes the implicit functions of the binary variables into a binary quadratic model using the factorization machine (FM) 28 and the Bayesian optimization of combinatorial structures (BOCS) 29 . This methodology is extended to applications such as the metamaterial design 24,25 , the chemical structure design 30 and the vehicle design 31 .
In the present paper, we present a QUBO formulation relying on BBO for structural design optimization. As an application, we deal with a printed circuit board (PCB) design in a vehicle for avoiding resonance. A PCB used in the power control unit of a hybrid vehicle is a key device to control electric power. The resonance of the PCB causes defects of the bonded electric parts and the connectors on the PCB. Therefore, the PCB tightens screws through the mounting holes so that the natural frequency of the PCB becomes high to avoid low-frequency resonance. One of the solutions to maximize the natural frequency is to fix the PCB with as many screws as possible. However, the number of mounting holes is restricted because there are many electric parts on the PCB, and the number of screws affects fabrication cost in mass production. Consequently, the PCB design requires maximizing natural frequency and minimizing the number of mounting holes. We formulated this design problem as a multi-objective optimization problem, which is solved by two classical methods, the weighted sum method and the ε-constraint method 32 .
We demonstrated the performance of the FM-based black-box optimization for a QUBO solver in the PCB design. Figure 1 illustrates the PCB design procedure. BBO approximates the natural frequency to a binary quadratic model followed by FMQA 24 . The natural frequency is a solution to the eigenvalue problem derived from the equation of motion. We call the calculation of the eigenvalue problem a frequency analysis. The eigenvalue problem and its solution were dealt with the constraint and objective in the optimization process. Then we employed FM as BBO to approximate the natural frequency obtained from the frequency analysis to a binary quadratic model. By doing this, the PCB design could be expressed in QUBO. In the present paper, the optimized solution to the QUBO and its neighbors which are probabilistically selected were utilized to update the BBO, while the original FMQA 24 did not explicitly utilize a solution other than the optimized solution. The QUBO constructed by FM (hereinafter, our algorithm based on FMQA is called FM-QUBO) was applied to design PCB models with 17 or 27 candidates of the mounting holes. As a result, the proposed FM-QUBO showed higher success ratio of finding an optimal solution and shorter calculation time than the QUBO constructed by BOCS (BOCS-QUBO) 27 . In addition, we confirmed that the selection of QUBO solvers, quantum annealing and simulated annealing, leaded to the difference of the success ratio.

Results
Multi-objective optimization for resonance avoidance in PCB design. The design of the mounting hole positions on the PCB to avoid resonance is formulated as a multi-objective optimization problem (MOOP). When the natural frequency f of the PCB takes a high value, this design objective is satisfied. In addition, the number of mounting holes N has to be minimized to save space on the PCB and reduce the production cost. Therefore, we formulate this design problem as a MOOP with two objectives.
The mounting hole positions were represented by binary design variables, x = {x 1 , . . . , x n } ∈ {0, 1} n . We prepared n candidate positions of the mounting holes on the PCB design model. Then, x i = 1 if position i is selected as the mounting hole, x i = 0 if not. The number of the mounting holes N can be calculated as the sum of x i :

Estimation Simulation Optimization
Approximation to binary quadratic model using FM/BOCS Frequency analysis based on FEM Exploration of optimal solution by QUBO solver www.nature.com/scientificreports/ The natural frequency f depends on the mounting hole positions, i.e., the design variables x . In the frequency analysis based on the finite element method, f is calculated as a solution to an eigenvalue problem: where is the eigenvalue that is related to the natural frequency, i.e., = 2πf 2 , u is the eigenvector, and K and M are the stiffness and mass matrix calculated from the finite element model of the PCB, respectively. In the present paper, we only dealt with the lowest natural frequency. During the optimization, we assigned boundary conditions for fixing the displacement of a PCB model to the mounting hole positions corresponding to x . Then, since K and M depend on where the boundary conditions are assigned, they have to be reassembled whenever x changes. Therefore, G(x) is solved together with the optimization. The MOOP for the PCB design can be formulated as We regard the eigenvalue problem, G, defined in Eq. (2), as a constraint in the optimization. A QUBO solver cannot solve this design problem directly because f and G is an implicit function of x . Then, for expressing f as the binary quadratic model of x explicitly, we employed a machine learning technique and a Bayesian treatment, the factorization machine (FM) 28 and the Bayesian optimization of combinatorial structure (BOCS) 29 . Both approaches estimated the coefficients and constants of the binary quadratic model, Q ∈ R n×n , b ∈ R , based on data sets of x and f, which are the inputs and outputs of the eigenvalue problem G. The approximated natural frequency f is defined as For further details of the approximation calculation, see "Methods" section . Using f , the optimization problem in Eq. (3) can be expressed by the quadratic polynomial of x without the constraint: We have to convert the MOOP in Eq. (5) into QUBO, namely, two objectives into a single objective function F(x): The goal of the MOOP is to search for a set of optimal solutions, called Pareto-optimal solutions, which have trade-off relations. We employed two classical multi-objective optimization methods: the weighted sum method and the ε-constraint method. These methods scalarize multi-objectives into a single objective.
In the weighted sum method, the sum of the objectives with a weight w ∈ [0, 1] is utilized as the objective function: where α, β are parameters to normalize the objectives, N and f . Equation (7) has to be solved with different w to obtain Pareto-optimal solutions. This method can easily scalarize a set of objectives. However, the drawback is that all optimal solutions cannot be found if the MOOP is non-convex 32 .
The ε-constraint method can be applied to the non-convex MOOP. In this method, the objective functions except for one objective function deal with a constraint function, that is, When the MOOP in Eq. (8) is solved with all possible values of the parameter N ∈ Z , all Pareto-optimal solutions are found. The above problem is converted into a QUBO with a penalty parameter p: where p takes an appropriate value so that the second term becomes larger than the first term. One can also select the natural frequency f as the constraint and minimize the number of mounting holes. However, we considered that the formulation in Eq. (8) would suit the PCB design problem because the constraint value N is an integer value, and all Pareto-optimal solutions are found in n times optimization.

Algorithm of FM-based black-box optimization for QUBO solver. The algorithm of FM-QUBO
in a PCB design proceeds as shown in Algorithm 1. The algorithm basically follows FMQA 24 . One prepares a PCB model for the frequency analysis with sufficient initial and boundary conditions and assigns , are prepared in advance to calculate the PCB model. Then, the binary quadratic model in Eq. (4) that relies on FM is constructed to approximate d . In the following, the QUBO in Eq. (6) is solved by a QUBO solver and the optimized solution x * can be obtained. We, in the present paper, prepared three solvers, simulated annealing (SA), quantum annealing (QA), and random search (RS). After that, the frequency analysis is performed with x * to calculate the natural frequency f * , and d is updated. Moreover, we create two solutions x † , x † † based on x * and added d as described below. This process is iterated until the data set size m reaches m . Finally, a solution in d which minimizes the augmented objective function F (Eq. (7) or (9)) is regarded as the best solution for this design problem. We call these sequential procedures FM-QUBO and similar procedures using BOCS instead of FM BOCS-QUBO, respectively. Note that, in BOCS-QUBO, the neighbors were not used in the data adding process, i.e., m ← m + 1 , because BOCS is a probabilistic algorithm, which estimates the distribution of quadratic models of random variables and then samples a model.
In FM-QUBO, the neighbors of the optimized solution x † , x † † are added to the optimized solution. The data set size |d| at the beginning of the optimization is small. When the optimized solution is only added, an optimized solution strongly depends on the initial data sets. Then, we created two neighbor solutions x † , x † † of the optimized solution x * and added three data sets to d , m ← m + 3 in Algorithm 1, every time d is updated. It is expected that the neighbors of x * take similar or smaller objective function values. We define the neighbors based on the Hamming distance. One or two variables are selected from the optimized solution, and the selected variables are inverted, i.e., 0 to 1 or 1 to 0. This algorithm contributes to the search focusing on the solution space near the optimized solution, and the search for a global optimum solution due to its randomness.
In the present paper, the optimization was performed until the data set size reached m = 1800 . We conducted the calculation 20 times with different initial values of x . The average objective function values of the best solution and their 95% confidence interval were evaluated. Note that the best solution represents the solution that takes the smallest objective function values in the obtained solutions. In addition, we focus on how many calculations can find the optimal solution within 20 times calculation. Then, the probability of finding the optimal solution is defined as the success ratio.
PCB design based on the weighted sum method. We prepared a simplified PCB model (Fig. 2) to demonstrate how the present algorithms work while designing the PCB model using the QUBO formulated by  7) were set to 1/450 and 1/8, respectively. The parameter w was set to 0.5. We computed the optimal solution of this problem setting using a brute-force search, as shown in Fig. 3. The objectives of the optimal solution, i.e., the number of mounting holes N * and the natural frequency f * , are 8 and 825Hz, respectively. The optimal solution is one of the Pareto-optimal solutions. Other Pareto-optimal solutions will be found when w is changed in the range [0, 1]. At the end of the optimization, the objective function values of the best solutions of the FM-and BOCS-QUBO were not identical as shown in Fig. 4, while FM-QUBO required a shorter calculation time than BOCS-QUBO.   www.nature.com/scientificreports/ The FM-and BOCS-QUBO explored a smaller objective function value than the random search. In the early stage of the optimization process, i.e., the data size m was small, the natural frequency and the number of the mounting holes showed a different trend in FM-and BOCS-QUBO (Fig. 4b,c). As for the objective function value (Fig. 4a), BOCS-QUBO was converged slowly compared with FM-QUBO. However, when m became large, these values converged to a similar value regardless of the approximation method. Regarding the calculation time except for preparing the binary quadratic model, BOCS-QUBO took about 1.7 times longer than FM-QUBO. FM-QUBO found the optimal solution using smaller data sets than BOCS-QUBO. Figure 5 shows the success ratio for different data sizes of m = 250, 500, 1000, 1800 , respectively. When m became large, FM-and BOCS-QUBO reached the optimal solution in high probability. FM-QUBO obtained the optimal solution many times in the small m. At the end of the optimization, i.e., the data size m reached m = 1800 , the success ratios were 80% in FM-QUBO and 55% in BOCS-QUBO. The success ratio of BOCS-QUBO linearly increased with m. Therefore, BOCS-QUBO may need more data sets to acquire the optimal solution in higher probability.
PCB design based on the ε-constraint method. The simplified PCB model (Fig. 2) was designed using the QUBO formulated by the ε-constraint method, as shown in Eq. (9). We only show the result for the case of the constraint with the number of the mounting holes N of 6. The feasible solutions satisfying the constraint are 17 C 6 = 12,376 . The parameter p in Eq. (9) was set to 1/450. The optimal solution is illustrated in Fig. 6.
During the optimization, the average natural frequencies of the optimized solution of FM-and BOCS-QUBO increased with the data size m and demonstrated a similar trend, while that of random search gradually increased with m (Fig. 7a). At the end of the optimization, FM-and BOCS-QUBO could find the optimal solution of over 80% (Fig. 7b). In this problem setting, the success ratio of FM-QUBO was going up with m after m = 1000 , but that of BOCS-QUBO was not the same. This result shows an opposite trend compared to the results given in Fig. 5.
We prepared another PCB model with 27 candidates of the mounting holes, i.e., n = 27 , and additional masses, as displayed in Fig. 8a. This design problem was also formulated as the QUBO based on the ε-constraint method with N = 8 . The number of combinations of the mounting holes is 2 27 = 134,217,728 , and the feasible solutions are 27 C 8 = 2,220,075 . In this large-scale problem, the difference in the best natural frequency between FM-and BOCS-QUBO could not be identified from the optimization history (Fig. 8b) as well as the result in the  www.nature.com/scientificreports/ small-scale problem (Fig. 7a). However, BOCS-QUBO took 35.9 times longer than FM-QUBO for preparing the binary quadratic model.

Discussion
We presented the FM-based black-box optimization for a QUBO solver to design the PCB for resonance avoidance. The PCB design was formulated as the optimization problem with two objectives, the maximization of the lowest natural frequency and the minimization of the number of the mounting holes. The multi-objective optimization problem was scalarized based on the classical multi-objective optimization techniques to formulate the PCB design as the QUBO: the weighted sum method and the ε-constraint method. In the scalarized optimization problem, a BBO method using factorization machine (FM) with the probabilistic data selection approximated the natural frequency which could not be directly expressed as the binary quadratic model. The QUBOs constructed by FM (FM-QUBO) and a conventional BBO method using Bayesian optimization of combinatorial structures (BOCS-QUBO) were applied to design the simplified PCB models with 17 or 27 candidates of the mounting holes. In our prepared design problems, FM-QUBO could find the optimal solutions in the high probability in the problems with both small and large solution spaces, and require shorter calculation time than BOCS-QUBO. The average performance between FM-and BOCS-QUBO (Figs. 4a, 7a, and 8b) were not identical, while the success ratio (Figs. 5 and 7b) and the calculation time were different. FM-QUBO explored the large solution space using the small data sets (Fig. 5), and the elapsed time to estimate an approximation function was shorter than BOCS-QUBO. In the case of a relatively small number of feasible solutions (Fig. 7b), FM-and BOCS-QUBO reached the optimal solution with a higher probability at the end of the optimization. The difference in the calculation time directly depended on the time required for the parameter estimation in FM and BOCS. The  To identify the effect of the proposed algorithm which is employed randomly selected neighbors of the optimized solution, we show the optimized results using a simple algorithm which only used the optimized solution when the data d was updated, i.e., m ← m + 1 , in Fig. 9 (FM, SA, pt=1). Under the same data size m = 1800 , the simple algorithm only found the larger objectives' solution than the proposed algorithm and can reach the optimum solution at once within 20 times calculation. Therefore, we confirmed that the proposed algorithm with the probabilistic data selection contributed to finding a minimum solution and increasing the success ratio.
Regarding the optimization problems with constraints, since the neighbors included the unfeasible solutions, in the worst case, one-third of the data sets for the approximation only satisfied the constraints, i.e., N = N . Nevertheless, the proposed algorithm can find an optimum solution with higher success ratio comparable to BOCS-QUBO. We consider that the adding solutions violated the constraint condition would help precisely approximate the objective function F including the constraint function in the QUBO formulation, and result in the efficient search for the optimum solution.
The selection of the QUBO solver will affect the optimized result. The average performances of SA and QA were not identical, but QA explored the optimal solution at a slightly higher probability than SA (Fig. 9). If this result only led to the QUBO solver's performance, the QA would be a suitable solver for black-box optimization for a QUBO solver (BBO-QUBO) with many design variables. However, since we observed the small difference in the success ratio between SA and QA in a few examples, we cannot conclude that QA is a better solver than SA because heuristic optimization methods, including SA, utilize randomness for their efficient search in nature. In addition, we dealt with the small number of the design variables in the demonstrated problems. Therefore, in our problem setting, it may be difficult to find the advantage of QA regarding the calculation time as shown in the reference 34 . When the design problem which has a large number of design variables is solved by the proposed method, the QA advantages will be discussed. www.nature.com/scientificreports/ Finally, our results also suggest that the BBO-QUBO assists in expanding the QUBO solver's engineering applications. The present paper employed the lowest natural frequency as the objective and was approximated to the quadratic model using BBO. There are other design objectives in the PCB design, such as avoiding a specified frequency, reducing the deformation, and so on. Then, the objectives are evaluated using lower eigenvalues and their eigenvectors. As the eigenvalues are solutions to the same eigenvalue problem, the presented BBO-QUBO algorithms applies to these designs by modifying the objectives and constraints to be approximated. In addition, by replacing the frequency anaysis with other simulations corresponding to the design problems, the BBO-QUBO would be widely used in engineering, although it is essential to confirm its performance in each problem.

Methods
Approximation of natural frequency for QUBO formulation. Factorization machine (FM). FM was developed for learning sparse data efficiently 28 . Let introduce a binary quadratic model for the estimation of the observed data y: where w 0 ∈ R and w = (w 1 , . . . , w n ) ∈ R n are the global bias and the strength of i-th variable. v i = (v i1 , . . . , v ik ) ∈ R k is the strength vector of interactions between i-th variable and the others. The model parameters w 0 , w, v in Eq. (10) are estimated using a given data set of an input variable x ∈ {0, 1} n and the corresponding output value y ∈ R . When there is enough data sets, the FM in Eq. (10) can approximate the interactions between the variables. In general, k is chosen sufficiently small because the number of the estimation parameters v becomes small and the computation time can be reduced, O(kn). However, our objective in using the FM is to construct a binary quadratic model from the data sets, namely, we are not supposed to set k to a small value.
Before solving the QUBOs formulated by the weighted sum method (Eq. (7)) and ε-constraint method (Eq. (9)) as shown in "Results" , we searched the proper value of k in the QUBOs. We applied the FM-QUBO with different k values, k = 9, 12, 15 , to the QUBOs and compared the average performances and the success ratios (Fig. 10). From the results, in our problem setting, we confirmed that k affected the optimized results and the best k value depended on the QUBO formulation. Therefore, we employed the best k in "Results", namely, k = 12 in the section "PCB design based on the weighted sum method" and k = 15 in the section "PCB design based on the ε-constraint method". Note that the k dependency was not observe in the original paper of FMQA 24 .
Bayesian optimization of combinatorial structures (BOCS). BOCS algorithm was proposed to solve the discrete optimization problem using scarce data. An acquisition function ŷ(x) , which is an approximation function of observed data y in BOCS, is defined as: The parameter α = (α i , α ij ) is estimated by Bayesian treatment (see the reference 29 for more details) using the data sets of x and y. As the Gibbs sampling is performed in the Bayesian treatment, the parameter estimation essentially consumes a long time, O(n 3 ) . In general, a sufficiently large sampling size is required to approximate observed data precisely, and the size depends on the number of the design variables, n. Frequency analysis. The frequency analysis was performed based on the finite element method to calculate the natural frequency of a PCB . From the equation of motion of a PCB discretized by the finite element method, the eigenvalue problem can be derived as shown in Eq. (2), i.e., K ∈ R n non ×n non and M ∈ R n non ×n non are the stiffness and mass matrix obtained when the equation of motion is discretized using finite elements. n non is the number of nodes of the discretized system. K and M consists of the material properties, including density, Young's modulus, and Poisson ratio. and u are the eigenvalue and the corresponding eigenvector. is related to the natural frequency f, i.e., = 2πf 2 .
For the frequency analysis, we employed the open-source finite-element analysis software, CalculiX 35 . The finite element model of the PCB, which we used in "Results", was discretized by a quadratic triangular shell element, and its total number of nodes and elements were 19,113 and 9,406, respectively. We only considered two materials corresponding to the substrate and the mass component in the PCB model. Domains, D fix i , D mass l , were placed on the discretized PCB model in advance to assign the mounting holes and the additional masses. In the optimization, the boundary conditions for fixing the displacement and rotation are assigned to D fix i where the corresponding design variable x i becomes 1. The domains D fix j , in which x j = 0 , is regarded as the substrate. The density in D mass l is set to a higher value of the substrate to represent the additional mass, while Young's modulus and Poisson ratio were set to the same value. Figure 11 illustrates the frequency analysis result of the optimal solutions, as shown in Figs. 3 and 6.
Both SA and QA are heuristic optimization methods, based on the analogy between the search process in the www.nature.com/scientificreports/ optimization and the physical phenomena. SA was performed on a classical computer, while QA on a quantum computer that is developed by D-Wave Systems Inc. 4 . We performed SA as implemented in dwave-neal 37 , and QA as in D-Wave Hybrid Solver 38 . In FM-QUBO and BOCS-QUBO, SA was performed 100 times and QA 3000 times in an optimization calculation. The best solution of the optimized solutions was regarded as the optimal solution and used for constructing the approximation model of FM or BOCS.